!	Routines and data storage for embedded atom model and target 

	module eamdata
	implicit none

	! Number of defined EAM atom types
	integer :: neam = 0
	! EAM potential form selector
	integer :: eam_form = 2
	! EAM parameters and type data
	integer, parameter :: NPARAMS(2) = (/ 7, 13 /)
	integer, allocatable :: eamflags(:,:)
	character*8, allocatable :: eamnames(:)
	real*8, allocatable :: eam(:,:), inputeam(:,:), eammass(:)
	! Cutoff constants - do we want to make these user-definable variables?
	real*8, parameter :: base_rc = 1.95d0, base_rn = 1.75d0      ! Will be multiplied by the equilibrium nearest-neighbour distance later on
	logical :: eam_cutoffs = .true.
	! Forces and energies - embedded and pair potential components
	integer :: local_natoms = 0
	real*8 :: eam_energy_phi, eam_energy_F, eam_energy_total
	real*8, allocatable :: eam_forces_phi(:,:), eam_forces_F(:,:), eam_forces_total(:,:)

	contains

	subroutine alloc_eam(newneam)
	implicit none
	integer :: newneam
	neam = newneam
	allocate(eam(neam,NPARAMS(eam_form)), inputeam(neam,NPARAMS(eam_form)), eamflags(neam,NPARAMS(eam_form)), eamnames(neam), eammass(neam))
	end subroutine alloc_eam

	subroutine eam_energy_and_forces(natoms, eamtypes, rij)
	implicit none
	integer, intent(in) :: natoms
	integer, intent(in) :: eamtypes(natoms)		! Integer array defining EAM atom types for each particle
	real*8, intent(in) :: rij(natoms,natoms,4)	! Vector and distance array for particles
	select case (eam_form)
	  case (1)
            call eam_baskes(natoms, eamtypes, rij)
	  case (2)
            call eam_mei(natoms, eamtypes, rij)
	  case default
            stop "Unrecognised EAM potential form."
	end select
	end subroutine eam_energy_and_forces

	! Potential Form 1 - Baskes et al.
	subroutine eam_baskes(natoms, eamtypes, rij)
	implicit none
	integer :: i, j
	integer, intent(in) :: natoms
	integer, intent(in) :: eamtypes(natoms)		! Integer array defining EAM atom types for each particle
	real*8, intent(in) :: rij(natoms,natoms,4)	! Vector and distance array for particles
	real*8 :: Sij(natoms,natoms)			! Screening function for each atom pair
	real*8 :: dSij_drij(natoms,natoms)		! Derivative of screening function for each atom pair
	real*8 :: rhobar(natoms)			! Total background electron density for each atom
	real*8 :: drhobar_drij(natoms,natoms)		! Derivative of background electron density for each atom pair
	real*8 :: rhoa(natoms,natoms)			! Atomic electron densities between atom pairs
	real*8 :: drhoa_drij(natoms,natoms)		! Derivative of atomic electron densities between atom pairs
	real*8 :: Fi(natoms), dFi_drhobar(natoms)	! Value of embedding function for each atom
	real*8 :: dFi_drhoa(natoms,natoms)		! Derivative of embedding function for each atom pair
	real*8 :: phi(natoms,natoms)			! Pair potential between atom pairs
	real*8 :: dphi_drij(natoms,natoms)		! Derivative of pair potential between atom pairs
	real*8 :: r0, E0, A0, alpha, beta, rhoscale, Z0 ! EAM parameters
	real*8 :: x, omxsq, omxcubed, phibar, dphibar, rc, rn, drhoa
	real*8 :: astar, phi_bar, Eu, dEu, pot, fvec(3)

	! Allocate force arrays (once only)
	if (natoms.ne.local_natoms) then
	  local_natoms = natoms
	  allocate(eam_forces_F(local_natoms,3), eam_forces_phi(local_natoms,3), eam_forces_total(local_natoms,3))
	end if

	! For now, assume that only one atom type is present in the system
	r0 = eam(1,1)  
	E0 = eam(1,2)
	A0 = eam(1,3)
	alpha = eam(1,4)
	beta = eam(1,5)
	rhoscale = eam(1,6)
	Z0 = eam(1,7)

	! Calculate cutoff function values for each atomic pair
	! [Mei et al, PRb, 43, 4653 (1991), eq 18a/18b]
	rc = base_rc * r0
	rn = base_rn * r0

	if (eam_cutoffs) then
	  do i = 1,natoms-1
	    do j = i+1,natoms

	      x = (rij(i,j,4)-rn)/(rc-rn)

	      if (rij(i,j,4).le.rn) then
		Sij(i,j) = 1.0d0
		dSij_drij(i,j) = 0.0d0
	      else if (rij(i,j,4).ge.rc) then
		Sij(i,j) = 0.0d0
		dSij_drij(i,j) = 0.0d0
	      else
		omxsq = (1.0d0-x)*(1.0d0-x)
		omxcubed = omxsq*(1.0d0-x)
		Sij(i,j) = omxcubed * (1.0d0 + 3.0d0*x + 6.0d0*x*x)
		dSij_drij(i,j) = omxcubed * (3.0d0 + 12.0d0*x) - 3.0d0*omxsq * (1.0d0 + 3.0d0*x + 6.0d0*x*x)
		dSij_drij(i,j) = dSij_drij(i,j) / (rc-rn)
              endif
	      Sij(j,i) = Sij(i,j)
	      dSij_drij(j,i) = dSij_drij(i,j)

	    end do
	  end do
	else
	  Sij = 1.0d0
	  dSij_drij = 0.0d0
	end if

	! Calculate background electron density for each atom, and atomic densities for each atom pair neighbour distance
	! Jelinek et al, Cond. Matt. arxiv:cond-mat/0610602v4
	!
	!  Eq 4 : Background electronic density around individual atom i
	!		     rho_i^(0)
	!	 rhobar(i) = --------- * G(Gamma_i)
	!		      rho_i^0
	!
	!	rho_i^(k) = SUM  rho_j^a(k)(rij) * Sij			(Eq 9a, for k = 0)
	!		    j<>i
	!
	!	  rho_i^0 = rho_i0 * Z_i0 * G(Gamma_i^ref)		(Eq 7)
	!
	!				[		(  rij	    ) ]
	! rho_i^a(k)(rij) = rho_i0 * exp[ -beta_i^(k) * ( ----- - 1 ) ]	(Eq 10)
	!				[		( r_i^0	    ) ]
	!
	!	    r_i^0 = nearest-neighbour distance in reference structure (EAM parameter 1, r0)
	!      beta_i^(k) = element-dependent parameters (EAM parameter 5, beta)
	!	   rho_i0 = element-dependent density scaling (should be an EAM parameter 6, rhoscale)
	!	     Z_i0 = First nearest-neighbour coordination number of reference system (EAM parameter 7, Z0)
	!	      Sij = cutoff function between atom, calculated above
	!	Neglecting higher-order terms (i.e. only k = 0 considered), Gamma_i = 0 and G(Gamma_i) = 1 (also for Gamma_i^ref)

	rhobar = 0.0d0
	drhobar_drij = 0.0d0
	do i=1,natoms-1
	  do j=i+1,natoms

	    rhoa(i,j) = rhoscale * dexp(-beta*((rij(i,j,4)/r0)-1.0d0))
	    rhoa(j,i) = rhoa(i,j)

	    drhoa_drij(i,j) = rhoscale * dexp( -beta*((rij(i,j,4)/r0)-1.0d0) ) * (-beta/r0)
	    drhoa_drij(j,i) = drhoa_drij(i,j)

	    rhobar(i) = rhobar(i) + rhoa(i,j) * Sij(i,j)
	    rhobar(j) = rhobar(j) + rhoa(j,i) * Sij(j,i)

	    drhobar_drij(i,j) = drhoa_drij(i,j)*Sij(i,j) + rhoa(i,j)*dSij_drij(i,j)
	    drhobar_drij(j,i) = drhobar_drij(i,j)

	  enddo
	enddo
	! If k = 0 only, then Eq 7 reduces to 'rho_i0 * Z_i0'
	rhobar(:) = rhobar(:) / (rhoscale * Z0)
	drhobar_drij(:,:) = drhobar_drij(:,:) / (rhoscale * Z0)

	! Determine embedding energy for each atom
	! Jelinek et al, Cond. Matt. arxiv:cond-mat/0610602v4
	!
	!  Eq 3 : Embedding energy per atom
	!
	!  F_i(rhobar(i)) = Ai * Ei * rhobar(i) * ln(rhobar(i))
	!
	!	        Ai = Energy scaling parameter (EAM parameter 3, A0)
	!		Ei = Sublimation energy (EAM parameter 2, E0)
	Fi = 0.0d0
	dFi_drhobar = 0.0d0
	dFi_drhoa = 0.0d0
	do i=1,natoms
	  if (rhobar(i).gt.tiny(0.0d0)) then
	    Fi(i) = A0 * E0 * rhobar(i) * dlog(rhobar(i))
 	    dFi_drhobar(i) = (A0*E0)*(1.0d0 + dlog(rhobar(i)))
	  else
	    ! Note: Derivative at rhobar(i) = 0 should be -INF
! 	    write(0,"(a,i4,a)") "WARNING - rhobar for atom ",i," is zero...."
	    dFi_drhobar(i) = -(A0*E0)
	  end if
	end do
	do i=1,natoms-1
	  do j=i+1,natoms
	    if (rhoa(i,j).gt.tiny(0.0d0)) then
	      dFi_drhoa(i,j) = A0 * E0 * (1.0d0 + dlog(rhoa(i,j))) * drhoa_drij(i,j)
	    else
	      dFi_drhoa(i,j) = -(A0 * E0)
	    end if
	    dFi_drhoa(j,i) = dFi_drhoa(i,j)
	  end do
	enddo

	! Calculate pair potential between atoms
	! Jelinek et al, Cond. Matt. arxiv:cond-mat/0610602v4
	!
	! Eq 12 : Pair potentials
	!
	!    phi_ij(rij) = phibar_ij(rij) * Sij
	!
	!		    1  [		      ( Zij		    )	   ( Zij		 ) ]
	! phibar_ij(rij) = --- [ 2 * E_ij^u(rij) - Fi ( --- * rho_j^a0(rij) ) - Fj ( --- * rho_j^a0(rij) ) ]  (Eq 13)
	!		   Zij (		      ( Zi		    )      ( Zj		         ) ]
	!
	!	  E_ij^u = -Eij * (1 + a_ij(rij)) * exp(-a_ij(rij))						  (Eq 14)
	!
	!			     ( rij      )
	!	   a_ij = alpha_ij * ( ---- - 1 )								  (Eq 15)
	!			     ( rij0     )
	!
	!	   E_ij = Sublimation energy (EAM parameter 2 for pure metals, E0)
	!      alpha_ij = Exponential decay factor for universal energy term (EAM parameter 4, alpha)
	!	   Z_ij = "Depends on the structure of the reference system" (EAM parameter 7, Z0)
	do i=1,natoms-1
	  do j=i+1,natoms
	
	    astar = alpha*((rij(i,j,4)/r0) - 1.0d0)

	    Eu = -E0 * (1.0d0 + astar) * dexp(-astar)

	    x = A0 * E0 * rhoa(i,j) * dlog(rhoa(i,j)) ! * (Zij / Zi)
	    phibar = (2.0d0 * Eu - x - x) / Z0
	
	    phi(i,j) = phibar * Sij(i,j)
	    phi(j,i) = phi(i,j)

	    ! Derivatives
	    dEu = E0 * astar * dexp(-astar) * (alpha / r0)
	    x = dFi_drhoa(i,j)   ! * (Zij / Zi)
	    dphibar = (2.0d0 * dEu - x - x) / Z0

	    dphi_drij(i,j) = dphibar * Sij(i,j) + phibar * dSij_drij(i,j)
	    dphi_drij(j,i) = dphi_drij(i,j)

	  enddo
	enddo

	! Calculate forces on atoms
	eam_forces_phi = 0.0d0
	eam_forces_F = 0.0d0
	do i=1,natoms-1
	  do j=i+1,natoms

	    ! From pair potential
	    fvec(1:3) = -dphi_drij(i,j) * rij(i,j,1:3) / rij(i,j,4)
	    eam_forces_phi(i,1:3) = eam_forces_phi(i,1:3) - fvec(1:3)
	    eam_forces_phi(j,1:3) = eam_forces_phi(j,1:3) + fvec(1:3)
	
	    ! From embedding function
	    fvec(1:3) = -(dFi_drhobar(i) * drhobar_drij(i,j) + dFi_drhobar(j) * drhobar_drij(i,j)) * rij(i,j,1:3) / rij(i,j,4)
	    eam_forces_F(i,1:3) = eam_forces_F(i,1:3) - fvec(1:3)
	    eam_forces_F(j,1:3) = eam_forces_F(j,1:3) + fvec(1:3)

	  end do
	end do

	eam_forces_total(:,:) = eam_forces_F(:,:) + eam_forces_phi(:,:)

	! Calculate total potential energy
	eam_energy_phi = 0.0d0
	eam_energy_F = sum(Fi)
	do i=1,natoms-1
	  do j=i+1,natoms
	    eam_energy_phi = eam_energy_phi + phi(i,j)
	  end do
	end do
	eam_energy_total = eam_energy_phi + eam_energy_F

	end subroutine eam_baskes

	! Potential Form 2 - Mei et al.
	subroutine eam_mei(natoms, eamtypes, rij)
	implicit none
	integer :: i, j, k, m
	integer, intent(in) :: natoms
	integer, intent(in) :: eamtypes(natoms)		! Integer array defining EAM atom types for each particle
	real*8, intent(in) :: rij(natoms,natoms,4)	! Vector and distance array for particles
	real*8 :: Sij(natoms,natoms)			! Screening function for each atom pair
	real*8 :: dSij_drij(natoms,natoms)		! Derivative of screening function for each atom pair
	real*8 :: rho(natoms)				! Total background electron density for each atom
	real*8 :: drho_drij(natoms,natoms)		! Derivative of background electron density for each atom pair
	real*8 :: rhoa(natoms,natoms)			! Atomic electron densities between atom pairs
	real*8 :: drhoa_drij(natoms,natoms)		! Derivative of atomic electron densities between atom pairs
	real*8 :: Fi(natoms), dFi_drho(natoms)		! Value of embedding function for each atom
	real*8 :: phi(natoms,natoms)			! Pair potential between atom pairs
	real*8 :: dphi_drij(natoms,natoms)		! Derivative of pair potential between atom pairs
	real*8 :: r0, Ec, phi0, alpha, beta, gamma, delta, c(0:5) ! EAM parameters
	real*8 :: x, omxsq, omxcubed, dphi, rc, rn, drhoa, rhoe
	real*8 :: fvec(3), term, adivb, rdivr, dpart1, dpart2
	real*8, parameter :: roots(3) = (/ 1.0d0, 1.414213562d0, 1.732050808d0 /), s(3) = (/ 12.0d0, 6.0d0, 24.0d0 /)

	! Allocate force arrays (once only)
	if (natoms.ne.local_natoms) then
	  local_natoms = natoms
	  allocate(eam_forces_F(local_natoms,3), eam_forces_phi(local_natoms,3), eam_forces_total(local_natoms,3))
	end if

	! For now, assume that only one atom type is present in the system
	! Parameters:  Ec(eV) phi0(eV) r0(Ang) alpha beta  gamma  delta    c0       c1      c2        c3        c4       c5

	Ec = eam(1,1)
	phi0 = eam(1,2)
	r0 = eam(1,3)
	alpha = eam(1,4)
	beta = eam(1,5)
	gamma = eam(1,6)
	delta = eam(1,7)
	c(0:5) = eam(1,8:13)
	rhoe = 1.0d0

	! Calculate cutoff function values for each atomic pair
	! [Mei et al, PRb, 43, 4653 (1991), eq 18a/18b]
	rc = base_rc * r0
	rn = base_rn * r0

	if (eam_cutoffs) then
	  do i = 1,natoms-1
	    do j = i+1,natoms

	      x = (rij(i,j,4)-rn)/(rc-rn)

	      if (rij(i,j,4).le.rn) then
		Sij(i,j) = 1.0d0
		dSij_drij(i,j) = 0.0d0
	      else if (rij(i,j,4).ge.rc) then
		Sij(i,j) = 0.0d0
		dSij_drij(i,j) = 0.0d0
	      else
		omxsq = (1.0d0-x)*(1.0d0-x)
		omxcubed = omxsq*(1.0d0-x)
		Sij(i,j) = omxcubed * (1.0d0 + 3.0d0*x + 6.0d0*x*x)
		dSij_drij(i,j) = omxcubed * (3.0d0 + 12.0d0*x) - 3.0d0*omxsq * (1.0d0 + 3.0d0*x + 6.0d0*x*x)
		dSij_drij(i,j) = dSij_drij(i,j) / (rc-rn)
              endif
	      Sij(j,i) = Sij(i,j)
	      dSij_drij(j,i) = dSij_drij(i,j)

	    end do
	  end do
	else
	  Sij = 1.0d0
	  dSij_drij = 0.0d0
	end if

	! Calculate background electron density for each atom, and atomic densities for each atom pair neighbour distance
	!
	!  Eqs 2 / 3 : Background electronic density around individual atom i
	!
	!	rho_i = SUM  f(rij)			(Eq 2)
	!		j<>i
	!			 l   c(l) ( r0 )l
	!	f(rij) = rhoe * SUM  ---- ( -- )	(Eq 3)
	!		 	0,5   12  ( r  )
	rho = 0.0d0
	drho_drij = 0.0d0
	rhoa = 0.0d0
	do i=1,natoms-1
	  do j=i+1,natoms

	    do k=0,5
	      ! SPEED Precalcualte (r0/rij) power terms
	      rhoa(i,j) = rhoa(i,j) + rhoe * (c(k) / 12.0d0) * (r0 / rij(i,j,4))**k
	      drhoa_drij(i,j) = drhoa_drij(i,j) - rhoe * (c(k) / 12.0d0) * k * (r0 / rij(i,j,4))**k / rij(i,j,4)
	    end do
	    rhoa(j,i) = rhoa(i,j)
	    drhoa_drij(j,i) = drhoa_drij(i,j)

	    rho(i) = rho(i) + rhoa(i,j) * Sij(i,j)
	    rho(j) = rho(j) + rhoa(j,i) * Sij(j,i)

	    drho_drij(i,j) = drhoa_drij(i,j)*Sij(i,j) + rhoa(i,j)*dSij_drij(i,j)
	    drho_drij(j,i) = drho_drij(i,j)

	  enddo
	enddo

	! Determine embedding energy for each atom
	!
	!  Eq 5 : Embedding energy per atom
	!		     [     alpha    ( rho  ) ] ( rho  ) alpha/beta
	!  F_i(rho(i)) = -Ec [ 1 - ----- ln ( ---- ) ] ( ---- )
	!		     [     beta     ( rhoe ) ] ( rhoe )
	!	
	!			      m
	!		+ 0.5 phi0 * SUM s(m) exp( -(sqrt(m)-1) * gamma )
	!			     1,3
	!
	!		  [ 				      delta    ( rho  ) ]
	!		* [ 1 + (sqrt(m)-1) * delta - sqrt(m) ----- ln ( ---- ) ]
	!		  [				      beta     ( rhoe ) ]
	!
	!		  ( rho  ) sqrt(m) * gamma / beta
	!		* ( ---- )
	!		  ( rhoe )
	Fi = 0.0d0
	dFi_drho = 0.0d0
	adivb = alpha / beta
	do i=1,natoms
	  if (rho(i).gt.tiny(0.0d0)) then
	    ! SPEED Precalc powers, log
	    rdivr = rho(i) / rhoe
	    Fi(i) = -Ec * (1.0d0 - adivb * dlog(rdivr)) * rdivr**adivb
	    dFi_drho(i) = (-Ec*(1.0d0-adivb*dlog(rdivr))*(adivb*(rdivr**(adivb-1.0d0)))) + (Ec*adivb*(rdivr**adivb)*(rhoe/rho(i)))
	    do m=1,3
	      term = 0.5d0 * phi0 * s(m) * dexp(-(roots(m)-1.0d0)*gamma)
	      Fi(i) = Fi(i) + term * (1.0d0 + (roots(m)-1.0d0)*delta - roots(m)*(delta/beta)*dlog(rdivr)) * rdivr**(roots(m)*gamma/beta)
	      dpart1 = (1.0d0+(roots(m)-1.0d0)*delta-roots(m)*(delta/beta)*dlog(rdivr)) * (roots(m)*gamma/beta)*rdivr**((roots(m)*gamma/beta)-1.0d0)
	      dpart2 = -(rdivr**(roots(m)*gamma/beta)) * (roots(m)*delta/beta)*rhoe/rho(i)
 	      dFi_drho(i) = dFi_drho(i) + term * (dpart1 + dpart2)
	    end do
	  else
	    Fi(i) = 0.0d0
! 	    write(0,"(a,i4,a)") "WARNING - rhobar for atom ",i," is zero...."
	    dFi_drho(i) = 0.0d0
	  end if
	end do

	! Calculate pair potential between atoms
	!
	! Eq 4
	!			 [	     ( r      ) ]     [        ( r	) ]
	!    phi_ij(rij) = -phi0 [ 1 + delta ( -- - 1 ) ] exp [ -gamma ( -- - 1 ) ]
	!			 [	     ( r0     ) ]     [        ( r0	) ]
	do i=1,natoms-1
	  do j=i+1,natoms
	
	    rdivr = rij(i,j,4)/r0

	    term = -phi0 * (1.0d0 + delta*(rdivr-1.0d0)) * dexp(-gamma*(rdivr-1.0d0))
	    phi(i,j) = term * Sij(i,j)
	    phi(j,i) = phi(i,j)

	    dphi = -phi0*(1.0d0+delta*(rdivr-1.0d0))*dexp(-gamma*(rdivr-1.0d0))*(-gamma/r0)-phi0*dexp(-gamma*(rdivr-1.0d0))*delta/r0
	    dphi_drij(i,j) = dphi * Sij(i,j) + term * dSij_drij(i,j)
	    dphi_drij(j,i) = dphi_drij(i,j)

	  enddo
	enddo

	! Calculate forces on atoms
	eam_forces_phi = 0.0d0
	eam_forces_F = 0.0d0
	do i=1,natoms-1
	  do j=i+1,natoms

	    ! From pair potential
	    fvec(1:3) = -dphi_drij(i,j) * rij(i,j,1:3) / rij(i,j,4)
	    eam_forces_phi(i,1:3) = eam_forces_phi(i,1:3) - fvec(1:3)
	    eam_forces_phi(j,1:3) = eam_forces_phi(j,1:3) + fvec(1:3)
	
	    ! From embedding function
	    fvec(1:3) = -(dFi_drho(i) * drho_drij(i,j) + dFi_drho(j) * drho_drij(i,j)) * rij(i,j,1:3) / rij(i,j,4)
	    eam_forces_F(i,1:3) = eam_forces_F(i,1:3) - fvec(1:3)
	    eam_forces_F(j,1:3) = eam_forces_F(j,1:3) + fvec(1:3)

	  end do
	end do

	eam_forces_total(:,:) = eam_forces_F(:,:) + eam_forces_phi(:,:)

	! Calculate total potential energy
	eam_energy_phi = 0.0d0
	eam_energy_F = sum(Fi)
	do i=1,natoms-1
	  do j=i+1,natoms
	    eam_energy_phi = eam_energy_phi + phi(i,j)
	  end do
	end do
	eam_energy_total = eam_energy_phi + eam_energy_F

	end subroutine eam_mei

	end module eamdata
